Screening effects in a density functional theory based description of molecular 
junctions in the Coulomb blockade regime 
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We recently introduced a method based on density functional theory (DFT) and non-equilibrium 
Green's function techniques (NEGF) for calculating the addition energies of single molecule nano- 
junctions in the Coulomb blockade (CB) regime. Here we apply this approach to benzene molecules 
lying parallel and at various distances from two aluminum fee (111) surfaces, and discuss the distance 
dependence in our calculations in terms of electrostatic screening effects. The addition energies 
near the surface are reduced by about a factor of two, which is comparable to previously reported 
calculations employing a computationally far more demanding quasi-particle description. 
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A key issue in the emerging field of molecular elec- 
tronics is the description of electron transport between 
nanoscale contacts, for which considerable progress has 
been recently achieved at the experimental level^. The- 
oretically, two limiting regimes can be distinguished, 
namely, coherent transport (CT) for strong coupling be- 
tween the molecule and the electrodes and Coulomb 
blockade (CB) for weak coupling. The CB regime is best 
described by stability diagrams, where frontiers between 
low- and high-conductivity domains in bias and gate volt- 
age coordinates are reflected by diamond-like shapes^ii^ii^. 
For a proper description of these diagrams, the energy 
difference between the ionization and affinity levels of 
the inserted quantum dot or single molecule (commonly 
referred to as addition energies Fiadd) has to be evaluated. 

The first-principle non-equilibrium Green's function 
(NEGF) methods^'^'^'^ combined with density functional 
theory (DFT), which have been successfully used for 
the CT regime, are not so straightforward to apply for 
electron transfer under CB conditions, since an integer 
charge is transferred and results in a relaxation of the 
electronic structure of the central molecule. In principle 
only a many-body approach provides a general solution 
to the latter proble m^d°d^ , and even quasiparticle cal- 
culations based on the GW approximation^^ were found 
to not fully capture the impact of local spin and charge 
fluctuations in the CB regime-^. The suitability of a 
standard DFT framework for electron transport in both 
the CB and CT regimes was also debatcd^"'-^'''. due to 
the self-interaction (SI) of electrons^ in a Kohn-Sham 
framework (KS) and the lack of a derivative discontinu- 
ity (DD)ii in the evolution of KS -eigenenergies. 

As outlined above a main source of discrepancy with a 
DFT description of weakly coupled nanostructures relies 
on the fact that the gap between the highest occupied 
(HOMO) and lowest unoccupied (LUMO) molecular or- 
bital eigenenergies in a single particle KS scheme does 
not match in general the total energy difference between 
the ground state and lowest charged states when the size 
of the HOMO-LUMO gap is finitei^. This mismatch has 
been recently addressed in realistic calculations of Eadd 





FIG. 1: (Color online) Geometry and shape of the applied 
gate potential for a benzene molecule lying parallel and 
weakly coupled to two Al fee (111) surfaces for a distance 
of 8 A. The profile of Vgate (taken from the differences in spa- 
tial resolution for calculations at lOV and OV) is shown with 
grey shading, with its maximum located in the black regions. 



with standard DFT techniques in three different ways: i) 
For metal particles of finite size, a modified KS gap has 
been introduced, where the energetic difference between 
the HOMO (or LUMO) for charged and uncharged clus- 
ters has been directly taken into account^^; ii) For the 
description of the HOMO-LUMO gap in Cgo-metal in- 
terfaces, the charging energy has been obtained by using 
a constrained DFT formalism^Q^ where the occupation 
of hand-picked orbitals can be defined as a constraint in 
the input^i; iii) Within a NEGF-DFT framework, F,add 
has been defined via threshold values of an external gate 
voltage Y gate determined via a midpoint integration rule 
from induced charge transfer between small molecules 
(H2 and benzene) and lithium wires22.. 

In this work, we adopt the approach introduced in 
Ref.^^ and apply it to calculate ^add for a benzene 
molecule lying parallel and at various distances from two 
aluminum fee (111) surfaces (see Fig. [Ij. We focus in 
particular on the dependence of Ea^rf on the distance 
between the central molecule and the two electrodes, 
and argue that our method correctly describes screen- 
ing effects that lead to a reduction of the molecular gap 
due to the rearrangement of the electronic structure at 
the Al surfaces. This is supported by two key findings: 
i) the distance dependence of 'Eiadd (with a correction 
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FIG. 2: (Color online)Evolution of the added/removed elec- 
trons on the benzene molecule Nadd (obtained by spatial inte- 
gration of the electron density) as a function of Vgate for three 
distances between the molecule and Al surfaces (dA;s=4, 6 
and 8 A ). The values of Eadd reported elsewhere in this article 
are taken from such plots. 

term accounting for the geometric capacitance^'^) is co- 
herent with a purely electrostatic model for screening 
based on image charges^!; ii) The magnitude of screen- 
ing is comparable to values obtained for a variety of sys- 
tems with GW-^, constrained DFT^ and a recently de- 
veloped quantum-chemical approach^, where electrodes 
have been treated as a classical shape-dependent contin- 
uum, allowing to reproduce accurately the experimental 
data reported in Ref.-2.. 

Fig. [T] displays the system on which we performed the 
NEGF-DFT calculations with the commercially available 
ATK software^. The scattering region contains three 
layers of a 3x3 unit cell of Al on each side of the benzene 
molecule and three additionial layers on each side of the 
left and right electrode regions, respectively, where a 3x3 
k-point grid has been used for the sampling in the trans- 
verse Brillouin plane. All atoms in the Al layers have 
been left in their truncated bulk positions for the experi- 
mental lattice constant of 4.05 A. A double-zeta polarized 
(DZP) and single-zeta (SZ) basis set have been used for 
the molecule and Al surfaces, respectively, and LDA is 
chosen for the exchange-correlation (XC) functional. The 
Keldysh formalism^ allows for a self-consistent solution 
for the electron density of the open system as a whole for 
every value of the external gate potential V gat<r^- The 
shape of the effective potential generated by Ygate is il- 
lustrated as gray shades in Fig. [TJ 

In Fig. [2] we illustrate how we determine the addition 
energy from NEGF-DFT calculations, E^^'^'^, for varying 
aluminum-benzene distances Aaib, 

^add = lo dNaddVgate{Nadd) — J^i dNaddVgate{Nadd) = 

VgateiNadd = +0.5e) - Vgate{Nadd = -0.5e). (1) 

The first line in Eq. [T] is general in electrostatics and 
makes use of the relationship between a total energy E, a 



voltage V, and a charge N. In the second line, explicit 
threshold voltages obtained from NEGF-DFT calcula- 
tions are introduced, and E[^^'^ is evaluated from the 
gate voltages required to add or subtract half an elec- 
tron to the molecule. The relevance of Nadd=±0.5 comes 
from the integration only; the midpoint rule applies ex- 
actly in view of the linearity of the dependence of Nadrf 
on Vgate in DFT-MiiliiS, and since we calculate E'^Jdd 
the energy corresponding to the integral of Vgate over 
a transferred charge of ±1. This energy represents the 
input required for the transfer of one electron from the 
molecule to the two electrodes or vice versa in terms of 
the external potential Vgate inducing this transfer. We 
stress that this method includes screening effects implic- 
itly, as evidenced by the fact that the results do depend 
on dAiB (see Fig.H]). 

One should keep in mind that E^^'j^ consists of a sum 
of two terms, the first related to the modified molecu- 
lar HOMO-LUMO gap in the junction, the second to 
the electrostatic capacitance of the metallic electrodea^^. 
The latter is usually referred to as the geometric capac- 
itance contribution to the charging energy in the litera- 
ture^^ and will be denoted Egeom in the following; it can 
be safely neglected in the analysis of CB experiments on 
single molecule junctions since it scales with d^zs/A and 
the area of the electrodes A is usually well above tens 
of /im^, whereas the distance dAiB between the molecule 
and the electrode surfaces is in the A range. This is not 
the case, however, in our calculations. Because we apply 
periodic boundary conditions to the electronic structure 
in the plane perpendicular to the transport direction, the 
finite charges are transferred from the molecule to the 
metal surface in each unit cell. This means that A is 
defined by only nine atoms in the plane and only dAiB 
is of similar size as in the experiments, and as a conse- 
quence Egeom reaches the same order of magnitude than 
the energetic contribution from the molecular gap. In 
order to make our results meaningful with respect to ex- 
perimental values, we define a corrected addition energy 
as 

rpcorr _ rpcalc rp _ rpcalc ^ 6 {AaIB ~ 3^o) /„n 
^add — ^add ~ ^geom — ^add ~ 2 2eoA 

where e denotes the elementary charge and ep the dielec- 
tric constant of vacuum. Since the system is equivalent 
to two capacitors in parallel, the charging energy Egeom 
contains an additional factor i . Our expression for Egeom 
is rather approximate in the sense that it amounts to re- 
placing the molecule by a third metallic electrode with 
the same area A as the source and drain electrodes. In 
Eq. [2]dy!i;B-xo accounts for the difference in position be- 
tween the planes of the nuclei and the electrons focal 
points due to spilling effects, which defines the position 
of the surface in any purely electrostatic (and therefore 
not atomistic) model and corresponds to the image plane 
in the model for screening we introduce below^S.. 

We collect in Table |T] the values for for three 

different distances Aaib and contrast them with calcula- 



TABLE L Corrected addition energies 'Ef°JJ for three dis- 
tances Aaib between the Al fee (111) surfaces and the central 
benzene molecule; Egeom has been subtracted for a mean- 
ingful comparison with experiments. The uncorrected values 
EadtF ^re also given in parentheses. We provide in two ad- 
ditional columns the corresponding data for systems with Al 
and Li atomic chains as electrodes in order to connect our 
discussion to ReL'^'^. The lattice constants have been chosen 
as a.Ai = 2.39 and aLi ~ 2.9 A and the areas of the unit cell 
perpendicular to the wires as Aai = 4aA!x4ayi! and ALi ~ 
3aLix3a_Li. The position of the capacitor planes used for de- 
termining Egeom is taken as xo=2.0 A for both Al fee and Al 
wire electrodes and xo=2.3 A for the Li wires in accordance 
with the differences in interlayer spacinga^S. All values for 
EZ7 and EZ'd' are given in eV. 




d-AlB 


Al fee (111) 


Al wire 


Li wire 


4 A 


7.27 (8.69) 


8.26 (9.25) 


8.78 (9.80) 


6 A 


8.68 (11.51) 


10.18 (12.16) 


9.82 (12.03) 


8 A 


10.02 (14.27) 


11.19 (14.16) 


10.99 (14.40) 



tions for the one-dimensional systems studied in Refi^. 
For the largest distance of 8 A, the results for Li and 
Al wires as electrodes come rather close to the limiting 
case (Eg^p = 11.54 eV, as calculated from total energy 
differences for charged and neutral benzene molecules^^) , 
whereas the presence of the surface induces a gap reduc- 
tion for Al fee. A decrease in d^/s reduces E'^°JJ for all 
three types of electrodes due to screening effects which 
are distance dependent; the HOMO-LUMO gap at 4 A 
represents 63 %, 72 % and 76 % of E°^p with the Al sur- 
faces, Al wires and Li wires, respectively. These numbers 
are comparable to those found for other system o^°i^^'^^ . 
Although screening is usually associated with the inter- 
action of charges with surfaces, a similar albeit quanti- 
tatively smaller effect can also be observed in Table H] for 
the wire electrodes. 

In order to validate the E^^'J {dAis) values provided by 
our approach, we have also estimated the corresponding 
numbers from an image charge model. For that purpose 
we define the molecular contribution to Eadd (the capaci- 
tative term Egeom does not enter the picture here, since it 
is a correction for the finiteness of the unit cell, whereas 
the image charge model assumes an infinite surface by 
definition) as 

Kir = im + 1) - EiN)) (EiN) E{N 1)) = 

E% + S(7V + 1) + S(iV - 1) - E% ~ E,„.een (3) 

with E(N±1)=E"(N±1)-HS(N±1). EO^p denotes the 
difference between the electron affinity and ionization 
potential of the free molecule and S the correction 
due to screening (i.e. the image energies associated to 
the charges that are calculated following the detailed 
recipe given in the supplementary information of Refi^). 
T,{N)=Q, since there are no screening effects when the 



FIG. 3: (Color online) Evolution of E'^°JJ as a function of 
dAiB, as extracted from NEGF-DFT calculations and com- 
pared to E^^^" (for its definition see Eq. |3}. E°ap = 11-54 
eV as calculated in Refi^. 

benzene molecule is neutral, because it does not exhibit 
any polar bonds. 

We compare in Fig. [3] the distance dependence of Eadd 
as obtained from NEGF-DFT calculations via Eqs. [1] 
and [2] to the results provided by the image charge model, 
i.e. we test the assumption 

E':2ridMB)^EZ7{dAiB). (4) 

The deviations between E^™^^*^ and E'^°JJ are less than 
20 percent of the total value of for the range of 

dAiB values under consideration. This is quite remark- 
able given the approximative nature of: i) Escreen (defin- 
ing the dAiB dependence of E^™^^"^ via Eq. since only 
the charge distribution inside the molecule has been de- 
scribed realistically by using MuUiken charges^!; and ii) 
Egeom (contributing to E^^^'' via Eq. [5]) for which the 
molecule and metal surfaces have been replaced by ca- 
pacitor planes without taking into account any details of 
their atomic structures. When (Iaib gets smaller (close to 
4 A) , it has also to be considered that wavefunction over- 
lap, which is not included in electrostatic models, starts 
to play a role so that the agreement between E^™^^'^ and 
^add expected to be better at large distances. 

Finally, we want to position our work in the context 
of the other recently proposed methods for the theoret- 
ical description of CB experiments with single molecule 
jimctions. The modified KS scheme of Ref.^^ is based 
on finite systems and therefore not directly suitable to 
study screening effects. Our method differs from those 
in Rcfs.^^ and^^ by its level of accuracy; although infe- 
rior to a full quasi-particle description^^ , which can treat 
only rather small systems, our approach is preferable to 
a semi-empirical method^^, where the predictive power 
is limited by the need to find suitable parameters. The 
technique in Refi^ based on constrained DFT^- is rather 
close to our approach in the sense that it also enforces 
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the occupation of molecular orbitals and calculates Fiadd 
from the energy required to uphold this charging. How- 
ever, while we apply an external gate voltage Vgate, let 
the electron density relax as a function of it and deter- 
mine Eadd from threshold values for Y gate , the occupation 
of the HOMO/LUMO are fixed manually in Ref.-^, thus 
introducing a certain amount of arbitrariness (at least for 
not so weakly coupled systems), and the key quantity is 
the gradient of orbital eigenenergy with its occupation. 
We stress that all these methods find a gap reduction 
due to screening effects by about a factor of two, in good 
agreement with our results. 

In summary, we have extended a recently introduced 
method for the calculation of addition energies ^add for 
single-molecule junctions in the CB regime^- by describ- 
ing in a more realistic way the electrode surfaces. This 
paves the way towards NEGF-DFT based predictions of 



Eadrf for junctions characterized in recent experimental 
studies'^, where screening effects are likely to play a ma- 
jor role. We analyzed the distance dependence of Eadrf 
in comparison to an image charge model and to other 
techniques developed for determining Eadd and found an 
overall good agreement between all approaches suggest- 
ing a reduction of the electronic gap by up to 50 % in 
molecular junctions. 
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